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1.0  Introduction 

The  primary  goal  of  the  processing,  data  reduction,  and  scientific 
analysis  of  solar  wind  data  obtained  from  the  M.I.T.  experiments  on 
SOLRAD  1 1 A and  1 1B , was  to  provide  real-time  solar  wind  parameters  to 
DOD  users.  The  secondary  goal  was  an  extraction  of  detailed  parameters 
describing  the  solar  wind  in  the  magnetosheath  as  well  as  in  the 
interplanetary  medium.  It  was  hoped  that  special  attention  could  be 
given  to  periods  of  special  scientific  interest  and  that  the  high  time 
resolution  mode  of  the  instrument  could  be  used  to  advance  our  knowl- 
edge of  the  nature  of  magnetohydrodynamic  processes  in  the  solar  wind. 
In  particular,  the  first  month  of  operation  seemed  to  hold  special 
promise  for  interesting  physics  since  the  two  SOLRAD  11  spacecraft  were 
close  together  as  they  orbited  Earth.  An  additional  hope  was  that  the 
high  time  resolution  data  could  be  used  in  combination  with  radio 
scintillation  measurements  to  find  the  source  of  fluctuations  in  the 
wind  which  caused  the  scintillation  phenomenon. 

Most  of  the  anticipated  outcomes  were  partially  achieved;  but  as 
will  be  discussed  below,  the  potential  of  the  experiment  could  not  be 
fully  exploited . 

2.0  Real-time  data 

Preparations  were  made  before  launch  for  a rapid  analysis  of  the 
real-time  data  and  calculation  of  the  velocity,  density,  temperature, 
and  flow  direction  of  the  interplanetary  wind.  These  parameters  were 
to  be  sent  in  close  to  real  time  to  NOAA  where  they  were  to  form  a 
portion  of  the  SELDADS  data  base.  The  effort  got  off  to  a very  slow 
start:  the  programs  prepared  by  M.I.T.  were  ready;  but  the  cut  in 

funding  at  NRL  had  reduced  their  staff  to  such  a point  that  operational 
matters  had  to  take  precedence  over  data  analysis.  Although  SOLRAD  was 
launched  on  March  14,  1976,  there  remained  several  months  before  the 
problems  at  NRL  were  cleared  up.  The  data  were  finally  sent  to  NOAA  in 
1977,  and  the  problems  at  NCAA  were  not  solved  until  late  in  1977*  In 
the  intervening  period,  the  data  system  had  begun  to  scramble  our  data 
by  rotating  the  sequence  of  bits  which  represented  a complete  spectrum. 
A good  portion  of  our  effort  at  M.I.T.  was  expended  in  decoding  the 


rotated  data  and  developing  an  analysis  system  which  could  retrieve  the 
data  in  most  cases.  That  program  was  sent  to  NOAA  where  it  began  to 
function  normally.  NOAA  itself  was  plagued  with  a severe  funding 

shortage  and  much  to  do.  SOLRAD  had  a low  priority  late  1977.  The 

interaction  with  NOAA  required  our  clarifying  the  data  sequence  from 
NHL  for  them.  At  long  last  the  whole  system  worked.  Meanwhile,  SOLRAD 
11A  had  failed  totally  in  an  eclipse. 

3.0  Detailed  analysis  of  data 

The  full  description  of  this  process  and  the  scientific  results 
are  contained  in  a S.E.  thesis  by  Michael  Herrera.  The  thesis  is 

included  in  this  report  as  Appendix  A.  Briefly,  he  describes  the 

detailed  analysis  process  and  presents  results  from  different  regions 
around  the  earth.  The  determination  of  the  parameters  of  the  dis- 
tribution function  in  the  magnetosheath  is  of  considerable  interest. 
The  solar  wind  interaction  is  of  considerable  interest.  The  solar  wind 
interaction  with  Venus  can  be  characterized  in  a similar  way.  (Shefer 
1979). 

M.O  Periods  of, .specified  scientific  interest 

SOLRAD  solar  wind  data  analyzed  at  M.I.T  were  available  more 
rapidly  than  those  from  other  spacecraft.  They  were  used  for  diverse 
purposes . 

In  August  of  1976,  the  third  section  of  the  Skylab  workshop  on 
coronal  holes  was  held  in  Boulder,  Colorado.  Good  correlation  had  been 
established  by  the  M.I.T.  and  American  Science  and  Engineering  groups 
between  high-speed  solar  wind  streams  and  equatorial  coronal  holes 
viewed  in  soft  x-rays.  Solar  minimum  appeared  to  be  in  progress  as 
judged  from  the  absence  of  optical  and  UV  coronal  features.  It  was 
SOLRAD  data  which  showed  the  continuing  presence  of  a high  speed  stream 
and , hence , showed  that  the  relationship  to  coronal  processes  was  more 
subtle  than  had  been  assumed. 

The  relationship  between  coronal  holes  and  the  solar  wind  was 
studied  in  more  detail  in  the  paper  by  (Nolte  it  al. . 1977)  in  which 
SOLRAD  and  IMP  solar  wind  observations  were  used  in  comparison  with 
x-ray  coronal  images  to  show  that  high  speed  streams  were  very  much  in 
evidence  despite  a very  quiet  corona.  The  authors  concluded  that 
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either  coronal  holes  were  not  the  sources  of  all  recurring  high-speed 
solar  wind  streams  or  that  there  was  a change  in  the  appearance  of 
regions  which  gave  rise  to  the  open  magnetic  field  line  configuration 
thought  to  promote  the  emergence  of  high-speed  streams.  The  same 
observations  were  reported  at  the  American  Geophysical  Union  Spring 
Meeting  (Sullivan  et  al . . 1977). 

The  SOLRAD  data  were  also  of  considerable  assistance  in  estab- 
lishing a cross  calibration  between  the  Voyager  1 and  2 solar  wind 
instruments,  which  had  a ten  per  cent  uncertainty  in  calibration,  and 
the  SOLRAD  instrument,  which  was  calibrated  to  'v  three  per  cent. 

5.0  High-time  resolution  and  close  spacecraft  data 

A great  disappointment  in  the  progress  of'1  this  work  is  that  we 
were  unable  to  exploit  these  two  features  of  the  mission.  The  reason 
is  a mixture  of  lack  of  personnel  and  data.  It  should  be  clear  that  a 
good  deal  of  our  effort  was  expended  developing  an  algorithm  to  correct 
the  data.  Furthermore,  the  amount  of  high-time  resolution  data  was 
quite  limited  since  we  did  not  wish  to  (nor  were  we  allowed)  to  preempt 
time  from  other  experimenters  for  our  special  data  mode. 

In  practice,  that  meant  that  the  high-time  resolution  data  are 
only  available  near  the  end  of  the  SOLRAD  1 1B  lifetime  when  all  the 
other  operational  problems  of  the  spacecraft  settled  down.  We  cer- 
tainly cannot  object — it  is  clear  that  the  limited  funds  available  at 
NRL  forced  an  heroic  effort  to  keep  the  spacecraft  alive  at  all. 
Nevertheless,  there  is  still  potential  for  good  physics  in  the  data 
that  are  available. 

The  data  from  the  two  SOLRAD  spacecraft  taken  when  they  were  near 
to  one  another  are  also  potentially  valuable,  but  far  from  complete  as 
shown  in  Table  1. 


lafale  1 


19^7  - 05  - 238  11B  S/C  timing  circuit  failure 

missing 

97  - 116  312  - 327  Experiment  on  11A  failed  - Day  116 

11B  Data  bad,  Day  133 
11A  S/C  failed,  June  1977 

missing 

236  - 255  424  - 439 

missing 


296  - 472  - 

1978  - 20  - 543 

30  - 59  552  - 575 

Only  four  days  of  close  data  had  been  received  by  early  1979, 
therefore,  we  had  no  opportunity  to  explore  the  potential  of  the  data 
though  they  are  available  in  principle. 
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ABSTRACT 

This  is  a research  report  on  analysis  of  Experiment  15  on  the  Solrad  11A 
And  Solrad  11B  satellites  being  done  at  the  MIT  Center  for  Space  Research. 
This  experiment  involves  plasma  detectors  that  measure  current  within  differ- 
ential energy  channels  for  positive  ions  and  electrons.  The  satellites 
transmit  data  from  the  solar  wind  and  magnetospheric  regions.  This  report 
will  concentrate  on  a preliminary  analysis  of  the  positive  ion  distribution 
functions  determined  from  the  Solrad  data.  A graphics  display  minicomputer 
and  separate  Fortran  programs  serve  to  determine  characteristic  parameters 
of  and  plot  the  distribution  functions. 
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CHAPTER  ONE 

INTRODUCTION  TO  SOLRAD 
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Our  earth's  magnetic  environment  offers  an  ideal  opportunity 
for  the  study  of  a collisionless  shock  in  a plasma.  With  a long  mean 
free  path  and  relatively  long  Debye  length,  of  about  ten  meters,  of  the 
plama  surrounding  our  earth's  magnetic  field,  a satellite  is  an  ideal 
probe  for  surveying  the  characteristics  of  this  plasma. 

The  following  diagram  illustrates  the  basic  magneto spheric  regions 
about  the  earth  and  its  interaction  with  the  solar  wind. 


f t 


(diagram  taken  from  Vasyliunas) 

. The  solar  wind  is  a stream  consisting  of  protons  and  alpha 
particles  flowing  from  the  sun.  It  reaches  the  earth  with  an  average 
velocity  of  around  400  km/s  with  a density  of  about  10/cc. 

The  magnetopause  is  the  boundary  off  the  earth's  magnetic  field. 

It  is  located  at  a distance  where  the  pressure  due  to  the  solar  wind 
is  balanced  by  the  magnetic  field  pressure.  The  density  inside  of  the 
magnetopause  boundary  is  much  less  than  that  found  outside  of  it  as 
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the  magnetic  field  is  generally  impermeable  to  the  plasma  flowing  past  it. 

Since  the  solar  wind  is  hypersonic  at  the  earth  (with  Alv£n  wave 
speed  on  the  order  of  50  km/s),  a detached  shock  wave  exists  upstream 
from  the  magnetopause.  This  marks  the  outer  boundary  of  the  magnetosheath. 

The  magnetosheath  is  qualitatively  very  different  from  the  solar 
wind.  A collisionless  mechanism  at  the  bowshock  raises  its  temperature 
and  density  to  values  considerably  higher  than  the  solar  wind.  The  bulk 
velocity  is  reduced  and  its  angle  shifted.  Also,  there  exists  the 
phenomenon  of  a high  energy  tail  in  the  magnetosheath  spectra  that  differ 
from  the  Maxwellian  nature  of  the  solar  wind  spectra. 
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Experiment  J?  1 5 aboard  Solrad  11A  and  Solrad  1 IB  detects  ions  and 
electrons  of  plasma  in  the  solar  wind  and  magnetospheric  regions. 

These  satellites  orbit  in  a plane  close  to  the  ecliptic  at  a geocentric 
distance  of  20  earth  radii.  They  spend  approximately  equal  time  in  the 
solar  wind  and  in  the  magnetosphere. 

The  two  satellites  were  launched  simultaneously  and  remained 
close  together  initially  but  were  later  separated  by  180°  so  that  one 
spacecraft  is  always  in  the  solar  wind.  This  also  enables  more  complete 
tracking  by  tracking  sites  in  Blossom  Point,  Virginia  and  in  the  Arcetri 
Observatory,  Italy,  the  latter  of  which  has  only  one  antenna. 

This  report  is  a record  of  ongoing  research  at  MIT's  Center  for 
Space  Research. 


CHAPTER  TWO 


INSTRUMENTATION  - THE  FARADAY  CUP 

• 

At  the  heart  of  MIT’s  plasma  experiment  is  a Faraday  cup.  It  is 
comprised  of  a modulated  grid  situated  between  an  aperture  and  a collector 
plate  in  a cylindrical  package  that  points  to  the  sun.  The  grid 
modulation  selects  the  differential  energy  window  from  which  the  average 
value  of  the  distribution  function  is  calculated.  This  value  is  a function 
of  the  amount  of  current  striking  the  collector  plate  and  the  width 
of  the  energy  window. 


the  faraday  cup 

The  voltage  on  the  modulated  grid  is  a 500  Hertz  square  wave  that 

alternates  between  two  voltages  referred  to  as  VHI  and  VLOW.  The  DC 

current  from  the  collector  is  caused  by  positive  ions  with  an  energy 

per  unit  charge  associated  with  V , the  z - velocity  component,  that 

z 

is  greater  than  the  instantaneously  impressed  voltage  on  the  grid. 

The  difference  in  absolute  current  from  the  collector  between  that 
with  the  high  grid  voltage  and  that  with  the  low  grid  voltage  indicates 
the  number  of  particles  in  the  energy  range  defined  by  these  voltages. 
This  current  difference  is  directly  measured  as  an  alternating  current 
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180°  out  of  phase  with  the  modulating  voltage. 

The  collector  is  divided  into  three  sections  of  equal  area.  Utilizing 
the  ratios  of  currents  in  these  sections, the  angle  of  incidence  of  the 
plasma  relative  to  the  axis  of  symmetry  of  the  cylinder  of  the  Faraday  cup 
can  be  calculated  to  within  2°. 

Circuitry  controls  the  AC  and  DC  components  of  grid  modulation 
separately.  'Electrometers 'are  AC  coupled  to  each  collector  section. 

Any  problem  with  the  DC  photoelectric  effect  is  sufficiently  eliminated 
with  the  AC  coupling. 

Current  data  are  transmitted  from  the  spacecraft  by  means  of  data 

numbers  that  correspond  logarithmically  to  the  inverse  of  the  current  as 

follows  ( for  DN  the  datanumber  and  CURR  the  current  ) . 

DN  - 255.  + ln(6.94  x 10~13/CURR)/ln(l .0367) 

This  is  accomplished  by  a logarithmic  Analog/Digital  converter. 

In  this  manner  currents  through  a range  of  four  orders  of  magnitude 

can  be  encoded.  Each  data  number  is  written  as  a hexidecimal  pair  with 

-9 

00  corresponding  to  6.806  x 10  amperes  and  FF  corresponding  to 

6.94  x 10“  ^ amperes.  Because  the  current/data  number  conversion  is  logarithmic 

the  limiting  error  due  to  discrete  datanumbers  is  constant  throughout  the  range. 

In  our  experiment  24  contiguous  energy  channels,  detecting  the 
positive  ion  component,  range  from  241  ev  to  4279  ev  for  protons  and 
from  482  ev  to  8559  ev  for  alpha  particles. 

One  large  energy  channel,  the  30th,  measures  the  total  current 
from  ions  striking  the  collector  within  these  energy  limits.  The  currents 
from  each  individual  sector  of  no.  30  are  registered  in  channels  3 1,32, and  33. 

A second  Faraday  cup  positioned  90°  from  the  sun  measures  the  electron 
component  of  the  solar  wind. 
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Current  density  for  a collection  of  particles  each  with  individual 

charge  q and  overall  number  density  n equals  nqv  when  all  particles  have  the 

same  velocity  v.  The  current  collected  at  the  collector  of  the  faraday. 

« 

cup  must  therefore  be  equal  to  the  average  value  of  (nq$) • A(v)rs  where  A(v) 
is  the  area  of  the  portion  of  the  plate  that  is  illuminated  by  the  plasma 
and  ft  is  the  normal  to  it.  Thus,  taking  into  account  the  transparency,!, 
of  the  grids,  the  expression  for  the  current  is 


I - Tnq  :fcv.  f (v)A(v)dv 
*!  Z 


which  is  integrated  over  the  appropriate  limits  for  v • f(v)  is  the  velocity 

z 

distribution  function. 

To  determine  an  expression  for  f(v)  we  assume  that  it  can  be  considered 
constant  within  each  energy  window,  and  we  consider  the  angle  of  incidence 
and  plasma  temperature  small  enough  so  that  A(v)  can  be  considered  constant 
( the  aperture  area  is  smaller  than  the  collector  area  - in  Chapter  V the 
error  resulting  from  making  such  an  approximation  is  considered  ).  Then 
the  expression  for  fn,  the  distribution  function  value  in  the  nth  channel 
is  given  by 


for  I the  current  in  the  nth  channel  and  the  limits  corresponding  to  the 
range  of  the  channel. 

We  see  that  the  integral  expression  in  the  denominator  is  proportional 
2 2 

to  (v2n)  • (vl  ) and  therefore  to  the  energy  width  of  the  channel,  A E„.  As 

I am  considering  only  relative  values  for  f , normalization  is  not  considered 


o 


and  the  following  expression  is  used  without  the  determination  of  k: 


L I 


* a n 

n (q/m)AE 

n 

Plots  of  the  Distribution  Function 

The  Techtronics  4051  computer  was  then  used  to  make  plots  of  spectra 

for  the  data  transferred  from  the  original  data  tape.  The  abscissa  was 

scaled  to  show  v of  protons  with  each  tic  mark  corresponding  to  100  km/s. 
z 

The  spectra  could  then  be  seen  on  linear  or  semilog  plots. 

The  '15  Spectra'  program  inputs  the  datanumbers  for  15  cycles  from 
a tape  file  in  the  data  cartridge.  From  each  it  calculates  , from  the 
above  formula,  the  value  of  the  distribution  function.  These  values  are 
placed  in  ah  array  F(I,J)  dimensioned  as  F(24,15)  where  the  I corresponds  to  the 
channel  number  minus  one  and  J to  the  index  number  that  is  used  to 
reference  each  individual  spectra. 

The  VIEWPORT  function  in  the  program  allows  .the  graph  to  be  placed 
in  different  portions  of  the  screen.  The  following  pages  show  distribution 
functions  in  semi-log  plotted  for  various  regions:  the  solar  wind,  bowshock, 
magnetosheath,  and  magnetopause  regions.  The  parabolic  appearance  of  the 
solar  wind  plots  is  consistent  with  the  Maxwellian  function  being  known  to 
provide  a good  fit  for  these  cases.  The  magnetosheath  examples  cannot  be 
closely  fitted  to  a Maxwellian  and  Chapter  V shows  that  a Kappa  distribution 
may  provide  a very  close  fit. 

These  graphs,  and  others  I compiled,  seem  to  hold  much  information 
regarding  the  changing  of  plasma  parameters  in  the  bowshock  and  magnetopause 
regions  and  more  investigations  will  be  done  of  these  crossings.  One  problem 
to  overcome  is  the  separation  of  proton  and  alpha-particle  components. 

Another  factor  making  the  analysis  more  complex,  is  the  existence  of  multiple 


i 
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'shock  crossings.  Temporal  fluctuations  can  be  noticed  in  certain  spectra 

resulting  from  a bowshock  crossing  occurring  while  the  current  measurements 

are  being  taken.  How  rapidly  the  change  is  from  the  magnetosheath  to  the 

solar  wind  portion  of  the  curve  can  give  an  indication  as  to  the  shock’s 

thickness  dfvidedby  its  speed.  From  the  intended  analysis  of  shock 

velocity,  significant  values  of  thickness  and  speed  might  be  obtained. 

The  dots  placed  below  the  distribution  curve  represent  the  lowest 

value  of  f that  can  be  plotted  for  each  channel.  They  each  correspond 

to  the  datanumber  FF  which  signifies  that  the  original  currents  are 
-13 

6.94  x 10  amps  or  below.  The  numbers  on  the  top  left  hand  corner  of 
the  pages. containing  15  plots  refer  to  the  year  and  day  number  of  the 
year.  The  numbers  found  in  each  graph  are  the  hours.  The 
axes  have  the  following  scale. 


v /m/q  in  divisions  of 
Z Jm, fqf  100  km/s 

The  proton  and  alpha  particle  spectra  are  superimposed  on  the  graph,  i 
A computer  summation  of  a typical  case  with  two  Maxwellians  on  a semi-log 
plot  is  shown  below.  The  alpha-particle  peak  is  located  at  an  abscissa 
/2  times  greater  than  that  of  the  proton  peak  (as  is  expected  for  both 
components  having  the  same  bulk  velocity). 
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Page  10  shows  a series  of  solar  wind  samples.  Their  shape  is  the 
expected  sum  of  two  Maxwellians,  with  the  alpha-particle  peak  at  /2 
times  the  abscissa  of  the  proton  peak.  Bulkspeed  is  around  600  km/s.  The 
relative  density  (from  the  ordinate  being  qf/m)  of  alpha  particles  to  protons 
* 6%  by  number  assuming  both  components  have  the  same  thermal  speed. 

Page  11  shows  a series  of  magnetosheath  samples.  Thermal  speeds 
are  increased  so  that  identification  of  each  component  is  made  difficult 
(this  is  further  explained  in  chapter  5).  Magnetosheath  particles,  as 
is  observed  here  have  a characteristically  smaller  z-component  of  bulk 
velocity  (from  both  the  reduction  of  bulk  velocity  and  the  geometry  of 
the  shock).  The  z-component  of  bulk  velocity,  bulkvz,  is  in  this  sample, 

350  km/s. 

On  page  12  is  found  an  excerpt  from  a series  of  multiple  bowshock 
crossings.  The  solar  wind  speed  is  seen  to  be  about  400  km/s.  The 
'turbulence'  that  is  found  in  these  magnetosheath  spectra  seem  to  indicate 
a certain  superposition  of  particles  with  the  solar  wind  distribution.  An 
example  of  the  temporal  fluctuations,  mentioned  on  page  7 , seems  to  be  in 
the  spectra  at  hour  2323  at  the  300  km/s  position.  The  high-energy 
particle  fluctuation  in  2325  seems  intriguing.  It  might  signify  a periodic 
temporal  fluctuation,  between  magnetosheath  and  solar  wind  states,  having 
a period  on  the  order  of  1 sec. 


Considerable  information  is  unfortunately  lost  within  the  112.5 
seconds  between  plasma  measurements  near  a bowshock  as  many  shock  crossings 


would  pass  unnoticed.  The  fast  mode  should  be  utilized  more  often  to 
enable  virtually  all  crossings  to  be  observed  in  any  pass. 

A series  of  magnetopause  crossings  can  be  seen  on  page  13  . A 
decrease  in  density,  where  the  true  spectrum  often  falls  under  the  FF 
curve,  is  most  evident. 


Projections  of  three-dimensional  graphs  showing  the  distribution 
function  vs.  time  are  found  on  pages  14  to  17  . This  provides  an  informative 
perspective  for  viewing  temporal  (and  possibly  spacial)  changes  within 
each  energy  interval.  The  density  variation  in  the  solar  wind  spectra 
on  page  14  may  Indicate  a regular  compressional  wave.  In  the  magnetosheath 
sample  on  the  following  page,  regular  variations  also  seem  apparent. 
Although  a bit  more  difficult  to  visualize,  the  'plateau'  on  page  16  shows 
the  magnetosheath  spectra  in  the  middle  of  a double  crossing  of  the  shock. 
Page  17  shows  a magnetopause  crossing  where  the  peaks  along  the  time 
dimension  signify  the  magnetosheath. 

Two  parallel  lines  in  the  t-v  plane  indicate  the  minimum  and  maximum 
velocities  that  Solrad  can  detect,  215  km/s  and  905  km/s.  The  range  of  " 
the  time  domain  is  listed  under  'Solrad  11A'  or  'Solrad  11B'. 
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typical  example  of  a series  of  solar  wind  spectra  not  found  near  the  bowsho'ck 


A series  of  magnetosheath  spectra  found  away  from  the  bowshock. 
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CHAPTER  THREE 


SOLRAD'S  'ANALYSIS'  PROGRAM 


The  purpose  of  the  'Analysis'  program  is  to  return  values  of  the 
plasma  parameters  of  z-component  of  bulkspeed  (SPEED),  density  (DEN), 
most  probable  thermal  speed  (THSP),  and  the  temperature  (TEMP),  and 

provide  the  values  of  angles,  including  THETA,  the  angle  of  incidence, 

» 

involved  in  the  analysis,  from  the  original  data  tape.  The  WIND 
subroutine  of  this  program  provides  the  derivation  of  these  four 
parameters  by  rough  numerical  integration  of  the  zeroeth,  first,  and 
second  moments  of  velocity.  These  integrations  are  done  over  a 
significant  portion  of  the  distribution  function,  where  the  limits  must 
be  determined  by  the  position  of  the  peak  channel  of  the  distribution. 

The  variable  F(J)  used  by  the  program  is  proportional  to  a 'differential 
density'  of  the  plasma.  It  is  defined  as  the  current  in  the  Jth  channel 

divided  by  the  average  velocity  detected  by  the  channel.  That  is, 

WTx  CURR(J)  CURR(J) 

( ' = V(J)  Jj(vl+v2) 

for  vl  and  v2  the  velocity  limits  of  the  channel.  This  equals  the  density 
of  only  those  particles  within  the  energy  limits  of  the  channel  times  the 
charge  per  particle  times  the  area  of  collector  plate  illumination. 

The  distribution  function's  average  value  in  the  nth  channel  equals 

k CURR(J)  ^ CURR(J)  a F(J) 
n (q/m)AEj  Js(v2-vl)  (v2-tvl)  AV^ 

as  is  expected  from  the  definition  of  distribution  function  (f=3n/3v) 

and  F(J)  for  AE  the  energy  width  and  AV  the  velocity  width. 

J J 

The  program  wrongly  seeks  the  channel  of  highest  F(J)f  which  varies 
as  the  distribution  function  only  when  the  velocity  width  is  constant. 
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For  the  solar  wind  samples,  consecutive  channels  near  the  proton 
peak  vary  to  a greater  extent  than  in  a magneto sheath  sample.  This  large 
variation  overcomes  any  effect  from  not  considering  the  distribution 
function  and  none  of  the  solar  wind  spectra  seem  to  be  affected  by  the 
error.  In  the  magnetosheath,  however,  with  greater  thermal  speeds,  the 

distribution  function  varies  less  between  adjacent  channels  so  that  this 

» 

error  often  causes  the  computed  bulk  velocity  to  be  greatly  shifted  to 
the  right,  to  a velocity  within  or  near  those  associated  with  the  12th 
or  13th  channel  where  the  energy  width  increases  and  therefore  the  velocity 
width  changes  the  greatest  amount  from  the  channel  one  below  it  in  energy. 
Lines  30  to  31  of  the  program  determine  the  peak  of  F(J). 

The  region  of  integration  is  taken  to  be  from  three  channels  below 
to  three  channels  above  the  maximum  flux  channel  ( channel  M ) if  these 
channels  exist  ( see  line  35  ) . In  the  required  integrals,  the  correct 
summation  expression  involves  the  current  flux.  Consider  the  density. 


/ F(v)dv  « £f  Av  and  therefore  / F(v)dv  « EF(J) . 
n n 


Thus  it  is  the  following  summations  that  are  calculated  in  lines  36  to  46. 


DEN  = C5  EF(J) 

SPEED  * EV (J)F(J) 
C5  EF(J) 


TQTP  and  THSP  are  calculated  from  DEN,  SPEED,  and  the  second 
moment  E V^(J)F(J). 

The  E represents  summation  over  the  seven  chosen  energy  windows. 


When  the  analysis  is  applied  to  the  magnetosheath,  two  other 
problems  arise.  First,  the  maximum  flux  channel  is  found  too  close 
in  energy  to  the  lower  energy  limit  of  the  apparatus  so  that  the  M-3 
channel  does  not  exist  ( as  a proton  channel  ) . This  is  because  the 
plasma's  z-component  of  bulk  velocity  is  considerably  decreased  ( by 
approximately  a factor  of  two  ),  from  the  geometry  of  the  bowshock 

» 

and  its  mechanism  of  collision,  from  its  previous  value  out  in  the 
solar  wind.  The  analysis  program  does  not  attempt  to  determine  any 
parameter  and  returns  a -100.  for  each  as  can  be  seen  in  the  data. 

The  second  problem  with  magnetosheath  spectra  is  that  there  is 
a greater  error  in  the  density  calculation  because  the  seven  energy  channels 
are  not  as  sufficient  as  for  the  solar  wind  samples. 

To  provide  a means  to  test  this  analysis,  the  Solrad  'Currents' 
program  was  developed  so  that  currents  expected  in  each  channel  could 
be  derived  from  the  plasma  parameters.  In  the  next  two  chapters  can  be  seen 
the  comparisons  between  the  observed  currents  and  those  expected  from  the 
determined  parameters. 


CHAPTER  FOUR 


SOLRAD'S  ’CURRENTS'  PROGRAM 


Description  of  'Currents'  program 

The  currents  program  simulates  the  datanumbers  that  correspond  to 
the  twenty-four  contiguous  positive-ion  energy  channels  and  one  large 
energy  channel  from  the  given  information  of  plasma  z-component  df  bulk 
velocity,  most  probable  thermal  speed,  and  number  density. 

The  current  from  each  channel  is  computed  by  a three-dimensional 
integration  over  velocity  space.  As  we  presently  do  not  consider 
predicting  the  separate  currents  from  each  of  the  three  sectors,  the 
symmetry  of  the  detector  lends  itself  to  a two-dimensional  consideration. 
However,  this  program  was  adapted  for  Solrad  and  it  was  originally  required 
to  consider  the  three  velocity  components  separately. 

In  addition  to  calculating  the  current  expected  from  the  direct 
effect  of  grid  modulation,  the  program  takes  into  account  a refraction 
effect.  This  is  caused  from  a greater  particle  refraction  with  the 
impressed  grid  modulation  at  VHX  volts  than  at  VLOW  volts  (for  VHI  and 
VLOW  the  upper  and  lower  limits  of  the  square  wave  voltage). 

Without  the  effect  of  refraction! 


voltage  . VHI 


voltage  • VLOW 
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As  the  grid  voltage  varies,  so  does  the  collector  current,  at 

180°  out  of  phase.  The  refraction  effect  subtructs  a greater  amount 

of  collector  plate  DC  current  at  grid  voltage  VHI  than  at  grid  voltage 

VLOW  when  the  plasma  image  of  the  aperture  does  not  fall  entirely  on  the 

collector.  This  is  because,  at  a non-zero  angle  of  incidence,  the 

center  of  the  aperture  image  is  displaced  by  a greater  amount  from 

» 

the  center  of  the  collector  at  VHI  and  the  plasma  would  illuminate  a 

smaller  part  of  the  collector,  as  shown  in  the  previous  diagram. 

Integrating  v^  between  the  limits  of  the  considered  channel 

calculates  only  the  current  not  produced  by  the  refraction  effect.  To 

take  it  into  consideration  v must  be  integrated  out  to  positive 

z 

infinity,  since  all  velocities  greater  than  the  threshold  velocity 
(for  entering  the  grids)  can  contribute  to  the  refraction  effect.  The 
necessary  integration,  performed  numerically  in  this  program  is  as 
follows  (using  the  notation  from  page  5 ) : 


ln  " Tqn((  v2f(v)A(v)dl * * * VxdVydvz)v.vHr(  Hi  V (7) A(^)dVxdVydVZ) V-VLOlP 


for  vln  and  v2^  velocity  limits  in  the  nth  channel,  since  the  first 

term  corresponds  to  the  absolute  current  at  VHI  and  the  second  to  the 
absolute  current  when  the  grid  voltage  is  VLOW.  The  difference  of  these 

two  terms  corresponds  to  the  AC  component  of  current  expected  (CURR(NCHAN) 
for  NCHAN  - n). 

f(v)  Is  taken  to  be  a Maxwellian  specified  by  the  input  parameters. 
A(v),  the  area  of  the  collector  that  is  illuminated  by  the  plasma,  is 
calculated  from  the  geometry  of  overlapping  circles.  A full  description 
of  the  program  can  be  read  in  Appendix  A. 
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The  'Currents'  program  yet  does  not  consider  the  effect  of  alpha 
particles.  If  the  analysis  program  were  to  be  modified  to  determine  the 
proportion  of  alpha  particles  in  a spectra  and  roughly  their  thermal 
speed,  the  consideration  of  alpha-particles  by  'Currents'  would  be  an 
immediate  consequence. 

Computational  results 

To  determine  how  closely  the  'Analysis'  program  and  'Currents' 
program  are  inverse  'functions’,  the  WIND  subroutine  was  made  into 
'Wind2'  so  as  to  accept  the  datanumbers  punched  on  cards  by  'Currents' 
and  analyze  the  produced  currents.  By  inputting  into 

'Currents'  a set  of  parameters,  any  deviation  from  these  in  the  output 
of  'Wind2'  is  then  noted. 

Results.:' showed  that  the  velocity  correlation  is  phenomenal,  to 
within  .8%  when  thermal  velocities  reasonable  for  the  solar  wind  are 
considered.*  A fairly  large  discrepancy  is  found  with  density  and  thermal 
speed.  The  density  discrepancy,  as  expected,  increases  as  thermal  speed 
‘is  increased . For  v^  = AOO^anf  thermal  speed  equal  to  50^m^insity  is 
off  by  only  about  2%.^  However,  a thermal  speed  of  100  kra/s  reduces  the 
density  by  almost  a factor  of  two  and  the  thermal  speed  is  off  by  -38%. 
Apparently  the  'Analysis'  program  is  not  integrating  over  a sufficiently 
large  range  of  the  distribution  function.  A further  analysis  of  this 
problem  is  necessary  with  regards  to  modifying  the  analysis  program  to 
enable  a more  accurate  analyzation  of  the  magnetosheath  spectra. 

* Velocity  is  off  by  as  much  as  12%  for  thermal  speed  ■ 100  km/sec. 
t Thermal  speed  is  off  by  about  8%. 


Graphic  Results 

For  comparison,  Che  Techtronlx  4051  graphics  system  was  used 
to  plot  the  distribution  function  computed  from  the  datanumbers  produced 
by  'Currents'  on  the  same  graph  as  that  derived  from  the  original  data. 

The  following  pages  show  the  original  data  in  dots  or  (+)  and 
the  expected  distribution  from  the  'Currents'  program  with  a solid 
curve.  The  discrepancies  talked  about  in  Chapter  Three  concerning 
the  magnetosheath  samples  are  evident  from  these  plots.  That  the 
'Analysis'  program  was  not  designed  to  properly  analyze  magnetosheath 
data  is  apparent.  When  the  lower  energy  distribution  values  do 
not  vary  rapidly  enough  the  curves  are  found  to  shift  to  center 
around  the  12th  or  13th  channel.  This  can  be  noted  from  the  linear 
plot  on  page  32.  From  this  it  is  apparent  that  the  error  is  large. 

The  solar  wind  spectra  match  very  well  in  all  parameters  with 
a slight  discrepancy  again  observed  in  density. 
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CHAPTER  FIVE 


GRAPHIC  KAPPA  AND  MAXWELLIAN  FITS 
The  Techtronix  graphics  display  can  also  be  used  to  make  an 
analysis  of  the  spectra  by  plotting  a theoretical  distribution 
function  of  v^, whose  parameters  are  determined  by  trial  and  error, 
that  fits  the  observed  function.  » 

Each  parameter  of  a Maxwellian  or  Kappa-  distribution  function 
can  be  changed  by  using  a 'user  definable  key'  on  the  console. 

The  Kappa  velocity  distribution  is  of  the  following  form: 


r(K+l)  ,,  . |v-vbulk|2x-(ic+l) 
*ro^j"  {1,J  ' ) 


This  distribution  function  is  thought  to  provide  a good  fit  for 
spectra  in  the  magnetosheath.  The  Kappa  distribution,  which  is  a 
statistical  sum  of  Maxwellians,  has  one  more  adjustable  parameter 
than  the  Maxwellian,  the  value  of  kappa.  The  following  page 
6hows  a Kappa  speed  distribution  (where  the  velocity  distribution 
function  is  multiplied  by  |v-vbulk|2)  for  cases  of  kappa  equals  one 
to  five  and  infinity.  A Kappa  distribution  for  kappa  equals  infinity 
is  a Maxwellian.  The  graph  on  the  left  shows  each  curve  relatively 
normalized  while  the  one  on  the  right  sets  each  peak  at  a constant  value. 
As  is  most  evident  from  the  righthand  plot,  the  Kappa  distribution  has 
a high  energy  tail, although  its  velocity  distribution  is  symmetric. 

Since  in  Solrad  we  can  only  deal  with  the  z-component  of  velocity, 
a speed  distribution  analysis  is  not  possible.  The  velocity  distribution 
must  be  reduced  to  a one-dimensional  distribution  by  integrating  over 

all  v and  v . The  result,  without  calculating  the  normalization  C is: 

x y 


z (1+  vz~  )* 

tew* 

For  ease  in  fitting  the  data,  the  maximum  is  constrained  to  one 
value  (determined  by  one  user's  key).  Thus,  what  is  plotted  is: 


f(v  ) = C(l+V/(1+  -\)K 

Z K KW. 


Density,  although  the  parameter  is  printed  next  to  the  graphs,  was 
not  used  as  an  adjustable  parameter  in  obtaining  the  Maxwellian  or  Kappa 

I 

fits. 

Results 

The  results  are  promising.  Magnetosheath  examples  are  shown  in 
pages  40  to  45  . It  is  interesting  to  note  that  spectra  taken  near 
each  other  in  time  each  have  very  closely  the  same  relation  to  the 
closest  curve  obtained  as  a fit. 

I I 

The  day  198  hour  205  spectrum  on  page  42  shows  a magnetosheath 
sample  taken  from  very  near  the  magnetopause. 

Page  43  shows  a case  that  seems  to  have  kappa  less  than  two. 

■ 

Because  of  the  alpha  particle  peak,  the  spectra  fits  were 
attempted  around  the  proton  peak  and  the  highest  energy  points.  This 
is  only  justified  when  the  sum  of  alpha  and  curves  converges  with 
the  proton  curve  at  high  energies.  To  determine  whether  or  not  this 
convergence  exists,  the  program  was  modified  so  that  alpha  particles 

I 1 

would  also  be  considered.  Page  44  , day  182  hour  1101  shows  each 
component  separately  together  with  their  sum.  It  is  seen  to  have 
the  convergence  necessary  although  a better  fit  could  be  obtained. 
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A closer  fit  for  the  alpha  particles,  however,  seems  to  indicate  (page  45) 
a smaller  thermal  speed  which  would  conserve  the  convergence.  However 
it  is*  seen  that  with  a larger  alpha-particle  thermal  velocity,  the 
convergence  does  not  exist.  A deeper  analysis  of  this  convergence 
problem  is  warranted. 


It  is  worth  remembering  at  this  stage  what  initial  assumptions 
were  made  for  the  distribution  function  to  be  calculated  in  the  plotting 
programs.  It  is  important  to  check  that  that  the  error  resulting  from 
approximations  in  the  derivation  of  the  formula  used  does  not  so  greatly 
affect  the  accuracy  of  each  point  that  the  curve  fitting  analysis  would 
not  be  valid. 

For  this  purpose  the  'Currerr'  program  was  devised.  It  is  the 
'Currents'  program,  modified  so  that  the  AREA  array  alternates  between 
describing  the  exact  area  response  function  and  an  ideal  flat  response 
function  (where  each  element  of  AREA  = ANORM) . The  current  values 
are  measured  at  each  channel  with  each  response  function  and  the 
percentage  deviation,  of  the  currents  calculated  with  the  exact-area 
response  function  from  the  currents  calculated  with  the  ideal  response 
function,  is  determined. 

These  values  would  then  give  an  Indication  as  to  the  error  resulting 
from  considering  the  area  constant  (i.e.  a flat  response)  in  the  plotting 
program.  The  bulkspced  was  considered  at  500  km/s  at  normal  incidence 
on  the  detector  and  density  at  5/cc.  Thermal  speed  was  varied  between  75 
and  250  km/s  to  observe  the  relationship  between  the  error  of  current  in 
a channel  and  the  quantity  w/vz« 

As  we  consider  calculations  for  higher  thermal  velocities,  the 


i 
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largest  current  deviations  are  found  in  the  low  energy  channels  where 
it  appears  that  the  expected  current  is  very  low  but  the  refraction 
effect  is  relatively  high  (since  it  is  due  to  particles  of  energy 
greater  than  that  representative  of  the  channel).  There  exists  a very 
low  deviation  in  the  high  energy  channels  perhaps  since  the  refraction 
effect  involves  much  fewer  particles. 

The  negative  deviation  in  low  energy  channels  for  low  thermal 
speeds  is  yet  to  be  understood. 

What  is  evident  from  the  results  of  this  program  is  that  the 
current  deviation  in  the  high  energy  channels  is  very  small,  as  for 
example  to  within  1%  for  channels  16  to  25  for  the  150  km/s  thermal 
speed  case.  This  is  important  to  note  when  considering  distribution 
fits  through  the  high  energy  points.  The' low  energy  points  may, 
owever,  be  significantly  be  affected  by  choice  of  response  function. 

The  'Currerr'  program,  it  should  be  kept  in  mind,  might  need 
to  be  revised  since  the  AREA  array  is  not  the  only  thing  in  the  program 
that  sets  a limit  to  the  response  function  actually  determined  in 
'Currents' . 
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APPENDIX  A 


SOLRAD  'CURRENTS'  PROGRAM 


The  beginning  lines  of  the  program  set  the  values  of  constants. 


DIMENSION  CURRI261 
!NTEliER«2  ON  (26) 
1NTE0ER*2  OAY.hOUR 
REAL* A NBULK.NTHERM 
DIMENSION  AREA (10001 
DIMENSION  VOLTLOI32) 
DATA  VOLTLO  /0. 0,240 
161. 5. 938. 2. 1 0 IS. 5.11 
2.283I.3.3072.8.33IA. 

DATA  VOLTHI  /0. 0.319 
139.7.1016.1,1171.7,1 
20,3067.3,3305.2,3548 
DIMENSION  XL  < 2a) ,WL ( 
DATA  XL/-0. 981560634 

1 ”0 . 367831 499 

2 ♦ 0 .58731 795A 

3 -0.981560634 

A -0.367831499 

5 *0. 58731795A 

DATA  WL/0.0A71 75336. 

1 0.233A92537. 

2 0.2031 67A27 , 

3 0.0471 75336 , 

A 0.233492537, 

5 0.203167427, 

DATA  XH/-3. 86972490. 

1 -0.947788391 

2 •1.59768264, 
DATA  WH/2.65855168F- 

1 2 . 6049231  OF- 

2 5. 1 6079856F - 
DIMENSION  VX<12>  «QT 1 
OATA  RC.RA  /2. 750.1. 


, V0:T‘U(32) 
.8.314. 7 , 398. 
71 . • < 1405.3.1 
4.  ;ri  . 4, 380  0 
.3.  i '..8,476. 
40',.  i . 1639.7, 
.3-  9.2,403 

24} . Full  2) ,WH 
.-0.  '-'.4117256 
.-  . 1 5233409 
, .0. '69902674 
,-0.  '04117256 
.-0.  25233409 
,♦0.  '69902674 
0.10-939326,0 
0.2-.  147046,0 
0. 140079329,0 
0.  K ‘>939326,0 
0.?'.  »147046,0 
0. 1'.'. 078329,0 
-3 . 02063703, - 
• -0.3l42i.0376 
♦2.27950706,. 
7.8.67368704F 
1,5. TO  1 35238E 
2,3. 9GS3905RF 
SO( it) ,XSQ( 12 
50/ 


7.473.6.549.9.630.0. 705.1.781.7.8 

640. 7. 1878.0. 21 15.5.2351 .7.2591 . 1 
.6,4045.0,240.5/ 

8,561.6,627.2,707.8,783.1 .859.5,9 

1876.1.2111 .2.2348.4.2587.0. 2826. 
3.1,4278.7.4279.5/ 

(12) 

,-0. 7 69902674, -0.587317954, 

,»0. 125233409. *0.36783 1499. 
,♦0.904 11 7256,* 0.98 1560634, 

.-0. 769902674,-0.587317954, 

,.0.1 25233409.. 0.36793 1499, 
,♦0.904 11 7256.* 0.98 1560634/ 

. 1600 78329, 0.203 167427, 
.249147046,0.233492537, 
.106939326,0.047175336, 

. 160078329.0.203)67427, 
.249147046.0.233492537, 

. 10b939J26, 0.0471 75336/ 
2.27950708.-1.59768264, 

.*0. 3 1 424  0 376. *0.94 778839 1 . 

3. 02U63703, *3. 88972490/ 

-5. J.90539058E-3.5. 16079856E-2. 

-1 ,3.701  J5236F.-1 .2.6049231  OF- 1 . 
-3.8.57368704E-5, 2. 658551 68E-7/ 

) 


DATA  C1,C2  /6.94E-13, 1.0367/ 


The  following  variables  are  inputs  to  the  'Currents’  program: 

DAY, HOUR  = time  spectrum  taken 
NBULK  = bulk  velocity  in  km/sec 

NTHERM  = the  most  probable  thermal  velocity  in  km/sec 
DEN  = the  density  of  the  plasma 
BETA  = shown  below  in  degrees 

ALFA  = shown  below  set  equal  to  zero  in  this  present  version 


original  a, 8 representation: 


\> 


far a day  cup 


a «■/,  subtended  by  v and  its  projection  onto  the  yz  plane 
0 *7,  subtended  by  -z  axis  and  projection  of  v onto  yz  plane 


CURR(26)  is  array  listing  current  values  in  channels  1 to  26  calculated 
by  the  program 

DN(26)  is  array  listing  datanumbers  calculated  after  each  CURR  array 
is  calculated 

AREA(lOOOl)  stores  the  function  relating  overlapping  area  of  solar  wind 

image  and  collector  plate  to  the  linear  displacement  of  their  centers 

VOLTLO(32)  , VOLTHK32)  VOLTLO (NCHAN ) and  VOLTHI  (NCHAN)  are  respectively 
lower  and  upper  voltage  limits  of  the  square  wave  grid  modulation  in 
the  NCHAN  .channel 


XL(24) ,WL(24) ,XH(12) ,WH(12)  are  sets  of  constants  used  in  the  numerical 
Integrations 

VX(1'  ) ,QTISQ( 12)  store  intermediate  values  in  the  integration 
RC-2.75  inches-  radius  of  collector  plate 
RA-1.5  inches-  aperture  radius 


1 IVERS  7132 

LOGICAL**  rW (5) 

CALL  HACK (TW) 

WRITE  (6.1020)  FW<l),Fw!n.IVE»S 

1020  FORMAT (69H1CURRENT  ON  TH.  FULL  COLLECTOR  PLATE  USING  AN. EXACT  AREA 
l CALCULATION.  .T111.2A  • 27H  SOLRAO.  ENERGY  CHANNELS  ,T1U» 
28HVERSI0N  IA) 

ICOUNT  = 0 
00  1 I = 1.29 
1 CURR(I)  = 0.0 


IVERS  = date  of  revision  of  program  = 7132  *■  day  132  in  1977 


HACK  subroutine  at  LNS  computation  facility  returns  the  date  the  job  is  run 


CURR  array  set  to  zero 

the  following  values  are  used  for  AREA(lOOOl)  determination: 


PI=3. 1A1S927 

ANORM=PI*RA*RA 

TRC=2.»RC 

TRA=2.*RA 

RCSO=RC*RC 

RASO=RA*RA 

0NO1=(RC-RA)**2 

8ND2=RC*RC-RA*RA 

BN03= (RC*RA) **2 

HPSRSQ=0.5*PI*(RCSO*RASQ> 

VBAR  = -110. 

0=1. 6021E- 1 9 

REM=I.672S2E-27 

AN*10ES*0 


ANORM  is  aperture  area 

VBAR  = -110.  is  suppressor  grid  voltage  in  volts 
-19 

Q=  1.602  x 10  is  elementary  charge  in  coulombs 
-27 

REM  “ 1.67252  x 10  is  proton  rest  mass  in  kg 
AN  «=  lO^xQ  = value  used  in  normalization 


Determination  of  function  represented  by  AREA(lOOOl) 

The  area  of  overlap,  A,  shaded  in  the  following  diagram,  is  determined 
as  a function  of  U,  the  displacement  between  the  centers  of  the  circles 
A ■ area  in  sector  ADC  + area  in  sector  BDC  - area  in  ADBC 
A - (f^)ir(RC)2  + (|^)tt(RA)2 
A - a(RC)2  + b(RA)2  - xU 


- 2 0sxU) 


x ■ (RC)  sin  a ■ (RA)  sin  b 


cos  a ■ 


2 2 2 
P +RC  -RA 

2U(RC) 


cos  b 


ARG1  * cos  a ■ (U2+BND2)/(TRC-U) 
ARG2  « cos  b « (U2-BND2)/(TRA-U) 
from  * above: 


2 2 2 
U +RA  -RC 

2U(RA) 


} from  defined  terms 


A ■ c°s'1<-iS-8cRa2)iic2  + - RA  .into. 

- BA 

Therefore:  2U*RC 

A - RC2cos-1(ARG1)  + RA2cos~1(ARG2)-U(RA)sin(cos~1(ARGl)) 
This  equation  is  used  in  the  following  DO  LOOP: 


VELP«0*2.*Q/REM 

c •••CALCULATE  AREA  overlap  AS  FUNCTION  OF  DISPLACEMENT  of  centers  of 
C COLLECTOR  AND  PROJECTION  OF  APERTURE  ONTO  COLLECTOR 
00  5 1*2.10000 

U=2.*RA*FL0AT(I-| )/10000.*RC-Ra 
USO=U*U 

AR61 * (US0*9ND2) / ( TRC*U) 

ARG2* (USO-RND2) / ( TRa»U) 

AREACI)=HPSRS0-RCS0*ARSIN(ARG1)-RAS0#ARSIN(aRG2>-u*RA*SQRT<  u.-arg 

1? ) • < 1 • ♦ AKG?) ) 

5 CONTINUE 

AREA  < 1 ) =AN0RM 

AREA  11 0001 ) =0.0 

SLOPt = 1 0000. /TRA 

81  = 10000.*  <RA-RC)/TRA«1 .5 

C ••••DIMENSIONS  OF  METERS*  KILOGRAMS.  AND  SECONDS  - EXCEPT  FOR  CUP 
C WHICH  IS  IN  INCHES 
C •••  DEN  IS  DENSITY  PER  CC 

alf=o 
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The  Index  I specifies  10001  values  of  U in  constant  increments 
from  (RC-RA)  to  (RC-RA)  and  the  area  is  calculated  for  each  I. 

This  enables  addressing  of  AREA  to  get  a sufficient  approximation 
to  the  area  of  overlap  at  specific  abscissas  in  the  integration. 


Determining  INDEX  from  U: 

• INDEX  = SLOPE-U  + Bl  for  SLOPE  and  B1  given 

, . 10000  (U+RA-RC)  . . .. 

since  we  let  INDEX  = 2^ +1.5 

The  value  of  1 becomes  1.5  so  that  one  obtains  the  closest  approximation 
to  the  fixed  variable  INDEX  from  the  floating  terms. 

ALFA=ALF*0 .01745329 

7 READ (5,1200, END=7  7 77 ) DAY , HOUR , NBULK , DEN , NTHERM , BE 
1200  FORMAT (13,15, F7 . 1 , 3F5 . 1) 

8 ICOUNT  = ICOUNT+1 

IF  (ICOUNT, EQ. 64. OR. ICOUNT. EQ. 128)  GO  TO  8 
C ***  NBULK=  RADIAL  VELOCITY=Z-COMPONENT  IN  KM/SEC 
THERM  = NTHERM*1000 
BETA  = BE 

B0=BETA*0. 01745329 
COSAL=COS (ALFA) 

BULKVZ='  NBULK*1000 


ICOUNT  ■ an  index  thatnumbers  the  spectra 

If  ICOUNT  - 64  or  128  the  values  are  not  used  since  these  values 
and  0 specify  the  calibration  modes 

At  the  end  of  inputting  all  cards  with  parameter  data  at  line  7 

control  is  transferred  to  line  7777  where  'NORMAL  EXIT'  is  written  and 

the  program  stops. 

THERM  ■=  NTHERM- 1000  “ thermal  speed  in  m/s 

BO  “ .0174  BETA  =6  in  radians 
ALFA  “.0174  ALF  = a in  radians 


BULKVZ  » NBULK* 1000  = z component  of  bulk  velocity  in  m/s 
BULKVY  = BULKVZ  * TAN(BO)  = y component  of  bulk  velocity  in  m/s 
BULKVZ  “ 0.  ■ x component  of  bulk  velocity 

The  'Currents'  program  is  adapted  and  was  originally  designed  for 
a three-dimensional  problem.  It  can  be  essentially  reduced  to  a two- 
dimensional  integration. By  setting  a to  zero,  our  problem  is  a 
two-dimensional  one  with  the  value  for  8 corresponding  to  the  value 
of  THETA  in  the  Solrad  ’Analysis'  Program  where  the  angle  is  computed. 

The  following  diagram  shows  the  velocity  components  for  non-zero  a: 


VX(I)  “ XH(I)  * THERM  + BULKVX  are  the  x-components  of  velocity 
used  in  tie  Hermite  numerical  integration  method  described  in  page  56. 
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CONST  = (6.45E-<O*0.735«AN*(ALPHa/PI>#*(1.5)*DEN 
C •••  VBOT  AND  VTOP  ARE  LIMITS  OF  INTEGRATION  AS  SET  BY  THE  DISTRIBUTION 
C FUNCTION 

VB0T=8ULKVZ-3.5*TMERM 
VTOP=BULKV7*3.5*THERM 
NCHAN  = 2 
CURMAX=0.0 
77  CONTINUE 

C • VLGRID  AND  VHGRID  ARE  LIMITS  OF  INTEGRATION  AS  SET  BY  THE  VOLTAGES  ON 
C THE  GRIDS 

VLOW= VOL TLO (NCHAN) 

VHl=VOLTHI (NCHAN) 

VLGRID  = SORT ( VELPRO*VLOR> 

VHGRIO  = SORT (VELPRO#VHI ) 

IF  (VLGRID. GE. VTOP)  GO  TO  170 
IL«=  1 

* IH=2* 

VZ1L  = AMAXl (VLGRID. VBOT) 

V71H  = AMIN! (VHGRIO, VTOP) 

IF (VBOT ,LE. VHGRID)  GO  TO  17 
IL  = 13 
VZ1H  = VBOT 

_ A -3/2 

CONST  = (6.A5  x 10  ) (.  735)AN(y)ir  J'*DEN  for  y - ALPHA 

is  used  immediately  after  the  triple  integration 

-A  2 

(6.A5  x 10  ) converts  sq.in.  used  in  AREA  to  m 

.735  = transparency,  the  fraction  allowed  through  grids  ^ 

AN  *»  10^Q  which  converts  density  units  from  ///cc  to  II /m 

-3/2  3/2  1/2 

yir  ' = (y/ir)  y » distribution  function  normali2ation  times 

constant  used  in  Hermite  integration 


DEN  *=  density  of  plasma  in  H/cc 

VBOT  = BULKVZ  - 3.5  * THERM  VTOP  » BULKVZ  - 3.5  * THERM 

The  v^  integration  involves  Gaussian  Integration  and  limits  are 
therefore  needed  for  the  distribution  function.  With  these  limits,  the 
fractional  deviation  from  the  exact  value  is  -7. A x 10  ^ , within  the 


accuracy  of  the  integration  method. 

VHI  and  VLOW  become  the  upper  and  lower  voltage  limits  in  the 
grid  modulation. 


VLGRID  = (VELPRO-VLOW)  * (2Q*VL0W/mp) 

VHGRID  = (VELPRO* VHI)  = (2Q-VHI/mp) 

These  are  the  velocity  limits  corresponding  to  the  voltage  limits. 


This  is  from  QV  » *S,npvz 


the  relation  between  V and  proton  velocity  limit. 
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If  VLGRID  < VTOP  at  line  170  SUMLO  « SUMHI  - 0.0  so  CURR(NCHAN)  - 0. 
in  output  since  too  few  protons  exist  in  distribution  function  above  VLGRID. 

The  a-integration  is  divided  into  two  parts  to  take  into  account  a 
second  effect  producing  a pulsated  component  of  current  (CURR(NCHAN) ) . 

Both  of  the  following  contribute  to  the  value  of  CURR(NCHAN) : 

(1)  The.  difference  in  number  of  protons  admitted  through  at  grid 

voltage  VHI  and  VLOW.  * 

(2)  The  difference  in  area  of  interception' of  plasma  image  with  the 

collector  plate  at  VHI  and  VLOW  due  to  the  refraction  by  the 
grid  potential  of  protons  entering  at  an  angle  to  the  cup's  axis. 

(1)  is  given  by  the  integral  of  the  distribution  function  between 
VLGRID  and  VHGRID  and  this  integration  is  indexed  by  1£I<12.  Integration 
of  (2)  is  indexed  by  13<I<24  and  involves  integrating  over  all  velocities 
above  VLGRID. 

VZIL  = maximum  of  (VLGRID, VBOT)  « lower  limit  for  integration  (1) 


That  is,  no  integration  is  done  below  VBOT  or  VLGRID.  Similarly  with 
VZIH  ■ minimum  of  (VHGRID, VTOP) , no  integration  is  done  above  VHGRID  and 
VTOP. 

If  VBOT>VGRID  we  see  that  IL  ■ 1 and  IH  ■ 24  meaning  that  I would 
vary  from  1 to  24  performing  both  integrations.  If  VBOT  VHGRID  integration 
(1)  is  unnecessary  due  to  a low  distribution  value  and  IL»13  and  IH=24 
so  that  integration  (2)  is  performed  with  the  lower  limit  equalling  VBOT. 
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Constants  used  in  Integration 


AVABL  and  AVABH  are  midpoints  in  the  region  of  interest,  AVBEL 
and  AVBEH  respectively  half  of  the  range  in  each.  They  therefore  have 
the  values  shown  in  the  program. 

To  accomplish  the  triple  integration,  there  is  a nest  of  three 


DO  loops,  the  outermost  loop  doing  the  integration,  begins  at  the 


DO  120  statement. 


17 


20 


22 

24 


AVBEL  = 0 . 5*  < VZ1 H-VZlL  > 
AVABL=0.S» <VZ1H*VZ1L> 
AVBEH=0.5*  JVTOP-VZ1H) 
AVABh=0.5*(VT0P*VZ1H> 

LOOP=0 

VALUE=0.0 

DO  120  I=IL.IH 

IF ( I .GE . 1 3) GO  TO  22 

VZ  = AVeEI.*Xl  (I)»AVABL 

AVBE=AVBEL 

GO  TO  24 

VZ=AVBEH*XL  < I ) * AVABH 
AVfc)E=AVREH 
VZ INV= 1 ,/VZ 
PART=0®2./(REM*VZ**2) 
VLPART=VL0W*PART 
IF  (L00P.E0. 1 ) VLPARTsVHI^ART 
IF(VLPART.GT.1.)G0  TO  120 
B = 2.55/(l.«S0RT(l.-VLPART)> 


1.22/<1.*SQRTU.-VBAR*PART)»»1.02 


B is  the  term  used  to  determine  displacement  due  to  refraction 


as  shown  on  page  56.  Variable  LOOP  has  value  0 or  1.  Integration  is 


done  once  for  each  value.  The  integrations  are  identical  except  that 
LOOP  • 0 determines  refraction  for  modulator  grid  voltage  ■ VLOW  and 
LOOP  = 1 determines  refraction  for  modulator  grid  voltage  = VHI. 
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2 

PART  » Q/(JsmpV2  ) VLPART  = VLOW-PART  (for  LOOP  = 0) 

* - VHI * PART  (for  LOOP  = 1) 

If  VLPART  > 1 DO  loop  continued  with  next  Index  value  since  for 
Q*VL0W  > or  Q*VHI  > JsmpV^  there  would  be  no  expected  current. 

Hermite  Integration 

This  is  used  for  the  vx  Integration.  • 

» * 12 
/ e-  f (x)dx'  - Ia;fef(xk) 
k-) 

where  xk  are  the  zeroes  of  the  Hermit e polynomial  for  n = 12.  are 
weight  factors.  The  array  elements  XH(K)  equal  xk  and  WH(K)  equals  u^. 
The  innermost  integral  is  of  the  form 

/ e~^  ^vx“vx ' ^ AREA (vx ' , vy ' , v ' ) dvx ' 

-oo  J 

lc 

for  vx  « BULKVX.  Therefore  we  let  x *=  ^vx  = (vx-vx' ) /THERM.  Then 

• 

dvx  »Y  dx,  this  constant  appearing  in  the  constant  term  in  page  52. 
xk  - XH(I)  = - (BULKVX-VX ( I ) ) / THERM  for  VX(1)  = vx'  = integration 
variable.  Finally, 

VX(I)  *=  XH(I)*THERK:+  BULKVX 

are  the  x-components  of  the  velocity  used  in  the  integration  as  in  page 


Gaussian  Integration 

This  is  used  for  vz  and  v^  integration. 


1 12 
/4f(x)dx  = Eu;lf(x^) 

-1  k-1 


I 


0 


9 


for  x.  the  zeroes  of  Legendre  polynomial  for  n=12.  We  want  to  evaluate 

v2 

functions  of  the  form  / f(v)dv  where 

vl 

v2  v2- (vl+v2)/2  (v2-vl)/2 

/ f(v)dv  = f f (v+(vl+v2)/2)dv  = S f (v+(vl+v2) /2)dv  . 

Vl  vl-(vl+v2)/2  -(v2-vl)/2 

v2 

Let  x *=  (v/({v 2-vl)/2).  Then  / f(v)dv  = 

vl 

1 12 

— 20~Vl  / f (x)dx  = 5j(v2-vl)  E a)  f(x  ). 

* -1  i*=l  1 1 


This  explains  the  derivation  of  integration  constants  on  page  54  . 

Calculation  of  refraction 

The  refraction  of  one  ion  as  it  passes  through  the  grids  causes  it  to 
traverse  a horizontal  distance  1 that  is  a function  of  tangential  velocity 


D - dj  + d2 


V « grid  voltage 


51 . 

1 m m d, 
P P 1 


d,  * v t,  + altl 
1 2 1 


Let  v_  «■  /(v2  + v2  ) = initial  tangential  proton  velocity.  Then 
i x y 


_v  ± /(v^) 

t,  ■ z z m 

1 E_ 

-QV/mpd1 


m d,v 

£q 1±/(1“  S ..a» 


Jjm  v2 
P 2 
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Because  of  time-reversal  symmetry,  for  t^  the  time  to  traverse 


t,  * mpd2Vz  (1±  /(l-QV/(5sm  v*)) . 
2 QV  p 


Therefore,  the  total  time,  t,  is 


t " tl  + t2  " ~qv~  d-Ai-QV/OmipVj)). 


The  sign  is  negative  so  that  t ■ D/v  for  V - 0.  Therefore, 

z 


1 - vTt 


mv  D(l-/(l-QV/!sm  v2) 
z p z 

QV  ( 1+/ (l-  QV/^smv2 ) 


_T  ( — ) 

ir  1+  (1-QV/Jjm  v2) 
z p z 


Let  B ■=  — + + 1.02 

l+/(l-QV/Jsm  v2)  l-/(  1-Q*  VBAR/ljm  v2  ) 

p z p z 

VT 

Then  1 ■ — B.  B is  shown  in  the  program  section  on  page  54. 
z 

HOT  = -V/» (RA*RC) /B 
T0P=-B0T 

AL 1M= AMAX 1 (BOT • BULK VY-2. p*TMERM) 

BLIM=AMIN1 (T0P,BULKVY»2.i*THERM> 

IF  (BUM.  LE.  ALIM)  GO  TO  1 BS 
AVBE1=0.5* (BLIM-ALIM) 

AVAB1  =0.5*  (BUM* ALIM) 

DO  25  L=  1 « 1 2 
0T1=VX <L> »VZINV 
0TI50(L)-QT1*0T1 

C •••  X = DISPLACEMENT  OF  SOLAR  WIND  IN  X 01 RECT ION  OUE  TO  REFRACTION 
C FROM  CHARGED  GRIDS 
X=8*0T1 
XSO(L)=X»X 
25  CONTINUE 

ZPART=ALPHA« (VZ-3ULKVZ) **2 
IF(ZPA«T.GE.70.)7PART=70. 

PHIZ=VZ*EXP(-ZPART) 

VALUE  1 = 0.0 

The  maximum  1 that  can  produce  a current  on  the  plate  is  1*>(RA+RC).  The 

velocity  corresponding  to  this  displacement  has  the  value  BOT  for  the  transverse 

—ROT 

component  (v  ) . Since  RA  + RC  = (-£ii£)B,  BCT  - -v  * (RA  + RC)/B. 

Vz  2 


v integration  limits 


Therefore,  v needn't  be  integrated  for  v < BOT  or  v < -BOT.  Let 

y y y 

TOP  *=  -BOT.  Thus  ALIM  and  BLIM  are  limits  of  integration  for  since 

the  value  of  the  distribution  function  is  too  low  outside  of 

(BULKVY-2 . 5THERM  <v  < BULKVY+2 . 5THERM) 

-4 

which  encloses  all  but  4 x 10  of  the  area. 

ALIM  = maximum  of  (BOT, BULKVY  - 2.5THERM) 

BLIM  = minimum  of  (TOP, BULKVY  + 2.5THERM)  • 

Let  A=BULKVY-2.5  THERM  and  B*BULKVY+2.5  THERM.  Possible  situations  are  depicted 
as  follows: 


(I)' 

• «§ 
} > 


flu*  ALIM 


m 


AUM- 


BLIM  Alim 


Only  cases  (5)  and  (6)  should  be  avoided,  which  is  accomplished  by 
IF(BLIM. LE. ALIM)  GO  TO  185  which  writes  the  values  of  ALIM  and  BLIM  and  sets 
CURR(NCHAN) =0. , and  proceeds  with  the  program. 


AVBEl  •»  . 5 (BLIM  - ALIM) 
AVAB1  - .5 (BLIM  + ALIM) 


P ! 

t 

I 


* 

f i 


i 


E 


These  constants  of  integration  are, as  explained  for  v , used  with 

z 


the  Gaussian  integration. 


The  DO  25  loop  determines  the  plasma  displacement  in  the  x-direction 


at  each  abscissa  (VX(L))  used  in  integration  and  places  its  square  in  array 
XSQ(L)  : 

XSQ(L)  *=  ( since  VZINV  - v -1  on  page  54.  . 

vz  * 

Function  of  v integrated 

PHIZ-,  e-T<V»in.KVZJ> 

z 

If  y(v2“BULKVZ)  2>70  it  is  set  equal  to  70.  This  is  to  avoid  a size 
error  in  exponentiation. 


c **•  VY .VZ  = VELOCITIES  IN  REFERENCE  FRAME  OF  SPACECRAFT 
DO  110  J=1.12 
VY=AVBE1*XL(J)«AVAB1 
0T2=VY»VZINV 
0T2S0=0T2*0T2 

C **•  Y = DISPLACEMENT  OF  SOLA*  WINO  IN  Y DIRECTION  DUE  TO  REFRACTION  FROM 
C CHARGED  GRIDS 
Y=B*UT2 
YSU=Y«Y 

YPARTr ALPHA* (VY-BULKVY) **2 
IF (YPART.GE.70.) YPART=70. 

PH1YZ=Ph!Z*EXP(-YPART) 

VALUE2=  0 • 0 
DO  100  K= 1 , 1 2 
DS0=XSQ(K) *YS0 

C •**  CHECK  TO  SEE  IF  TOTAL  DISPLACEMENT  IS  COMPLETELY  ON  OR  OFF  COLLECTOR  PLATE 
IF (DSO.GE.RND3) GO  TO  100 
IF(DSQ.lE.unD1)G0  TO  AO 
C •**  DISPLACEMENT  DUE  TO  REFRACTION 
OIPSE=SQRT (DSO) 

INDEX=Sl0PE*DIPSE»B1 


' 


i 

I 


! 


The  DO  110  statement  starts  the  second  loop  that  integrates  v^. 

The  (y-displacement  of  the  plasma  image)  squared  *=  YSQ  for  YSQ  = (v  /v  )2B2 

y z 

obtained  in  the  same  manner  as  XSQ(L). 

The  function  integrated  is  obviously  a function  of  v^.  It  is 


PHIYZ  - 


-Y(v  -BULKVZ) 2 -y(v  - 
e z e 1 y 


BULKVY) 


2 
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PHIYZ  = PHIZ  e ' ' y 


-Y(v  - BULKVY ) 2 


PHIZ  e 


-YPART 


and  numerical  integration  on  is  performed  for  each  value  of  PHIZ.  As 
before,  the  maximum  value  of  YPART  is  70. 

The  v^  integration  starts  at  the  DO  100  statement  each  time  the  loop 
is  entered  with  a different  value  for  YSQ.  Each  value  of  XSQ  is  added  to  it 
with  each  index  value  of  the  loop. 

DSQ  = XSQ(K)  + YSQ 

Thus,  the  square  of  the  total  displacement  is  known  for  each  value  of  v^. 

The  following  cases  for  DSQ  are  considered: 

I)  DSQ  j>(RC  + RA)2 
i.e.  DSQ  1 BND3 

Loop  is  continued  for  next  value  of 
k since  no  current  contribution  to  integral 

II)  DSQ  <_  (RC  - RA)2 
i.e.  DSQ  <_  BNDl 

INDEX  = 1 for. area  determination 

III)  (RC  - RA)2  < DSQ  < (RC  + RA)2 
The  AREA  array  is  used 


DIPSE  *=  /(DSQ).  Thus,  since  INDEX  = SLOPE  • DIPSE  + B1  at  each  value  of  DIPSE 
to  reference  AREA  value  at  each  displacement. 


The  following  section  completes  the  triple  integration. 


© 


r 


40 

90 

100 

110 

120 


GO  TO  90 
I NOE Xs 1 

VALUE2  a VALUE2  ♦ WM (K > • A«EA < INOEX) 
CONTINUE 

VALUE  1 a VALUE  1 • WL  < J)  ••’HI  TZ*V*LUE2 
CONTINUE 

VALUt=VALUE*WL  < I )*AVBP1*AV8E*VALUE1 
CONTINUE 

VALUE=CONST# VALUE 
LOOP=LOOP*l 


Statement  90  completes  Integration  since,  as  shown  on  page  55, 


/VY(  vx-V)  AREA(v' ,v' ,v' ) = E U).  AREA(VX(K),VY(J),VZ(D) 

00  x y 2 av  ^ k 

* K«* 

for  - WH(K)  and  AREA  function  determined  by  INDEX.  VALUE2  originally 
set  at  zero. 

VALUE 1 = VALUE  + WL(J)  • PHIYZ  • VALUE2 

o 


IX  IX 

- I ( Z WH(K)AREA(INDEX) ) 

*i 


-Y(VY(J)-BULKVY)2-fy  -v')2 

6 0 2»  Z V 


since  VALUE  ■ 0.  originally. 

VALUE  - VALUE  + WL(I)  • AVBE1  • AVBE  * VALUE1 


- Til  (Z  WH(K)AREA(VX(K),VY(J),VZ(I))))e'Y(VY(I)"BULKVY)) 

<•  JL  •'"<  - 

e-Y(VZ(I)-BULKVZ)2z 

which  is  the  required  equation. 

To  normalize,  correct  it  dimensionally,  and  take  into  account 
transparency,  and  elementary  charge: 

\ 

VALUE  - CONST  • VALUE 

for  CONST  - (6.45  x 10  *)  (.735)  (10®)Qy(ytt)“^^  as  described  in  page  52. 
AVBE1  and  AVBE  are  needed  as  shown  in  page  59. 


|i 
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I 


•Jtm 


o 

GO  TO ( 1^5 • 1 30 ) .LOOP 
125  SUMLO=VALUF 

IT  (VHGRIO.GE.VTOP>  GO  TO  ISO 

1L=  1 3 

IH=24 

. GO  -TO  20 

130  SUMHl=VALUE 
1 AO  CURR (NCHAN)  = SUMLO-SUMHl 

IF(CURR(NCHAN) .LT.CURMAX)  GO  TO  47 
CURmAX=CURR (NCHAN) 

47  CONTINUE 

IF  (NCMAN.FQ.26)  GO  TO  50 

IF (CURR (NCHAN) .LT.CURMAX*. 0001)  GO  TO  49  * 

48  NCHAN  = NCHAN  ♦ 1 
GO  TO  77 

49  MN=NCHAN»1 

00  51  I J=MN* 25 
51  CURR ( I J) =0.0 
NCHAN=26 
GO  TO  77 

! 

LOOP  = LOOP  + 1 changes  LOOP  ° 0 to  LOOP  = 1 or  LOOP  * 1 to  LOOP  ■ 2 . 


Then,  if  LOOP  = 1 goes  to  125.  For  LOOP  = 2 goes  to  130  since  both 


integrations  are  completed. 


LOOP  = 1 


Integration  (1)  and  (2)  ( page  53  ) just  performed  with  VLPART 


VLOW* PART  from  which  the  refraction  is  calculated  for  VLOW. 

Thusly,  integration  (1)  is  finished  since  no  particles  enter  grids  with  w 
velocity  less  than  VZIH  when  grid  voltage  = VHI. 


Integration  (2)  must  be  performed  again  to  subtract  DC  current  at 
voltage  = VHI  from  DC  current  at  VLOW  to  determine  the  modulation  resulting 
from  refraction. 


If  VHGRID  VT0P  there  is  no  need  to  perform  integration  this  second 

time  since  too  few  protons  will  enter  with  velocity  greater  than  VZIH. 

For  second  integration  IL=13  and  IH=  24  are  lower  and  upper  limits 

of  the  index.  After  these  are  set  program  goes  to  line  20  which  starts 

v integration  DO  loop, 

z 


SUML0  ■=  VALUE  for  LOOP  = 1. 


LOOP  - 2 


r 


Integration  performed  with  VLPART  - VHI-PART.  Then  finally 
CURR(NCHAN)  **  SUMLO  - SUMHI  since  we  wish  to  subtract  the  results 
from  each  integration. 


I 


' • 


CURMAX  has  the  value  of  the  CURR(NCHAN)  which  is  greatest  for  NCHAN 
from  2 to  its  present  value.  * 

If  NCHAN  •=  26  then  it  goes  to  50  where  the/ last  step  of  the  program 
is  to  calculate  the  26  data  numbers.  If  not  it  goes  to  77  and  sets  up 
the  next  integration. 

If  CURR(NCHAN)  < .0001  CURMAX  set  all  remaining  CURR (NCHAN)  to  zero 
since  all  these  currents  will  eventually  correspond  to  the  datanumber  FF. 

Control  is  transferred  to  statement  50  after  all  26  CURR(NCHAN)  are 
calculated. 

The  DO  57  loop  determines  the  datanumbers,  DN,  for  NCHAN  = 2 to  26, 
the  24  contiguous  positive  ion  channels  and  one  large  channel). 

The  DN  array  is  printed  first  in  decimal  notation  and  secondly  in 
hexadecimal  notation.  The  input  parameters  are  printed.  Day,  hour,  and 
DN  are  punched  on  cards  to  enable  the  WIND2  program  to  make  a check  on  this 
program. 


! 
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C ***  CALCULATION  FOLLOWS  FOR  DATA  NUMBERS 
50  DN(l)=ICOUNT 
DU  57  1=2,26 

IF(CURR(I)  .LE.6.94E-13)  GO  TO  54 
DN (I) =255 . 0+ALOG*C1/CURR (I ) ) /ALOG (C2) 

GO  TO  57 
54  DN(I)=255 
57  CONTINUE 

WRITE  (6, 1070)  DN 
WRITE (6, 1007) DN 

1070  FORMAT (1HO, 9X, 13, 2X, 24(13, IX) , IX, 13) 

1007  FORMAT (1H  , 10X, 22, 3X, 24 (Z2, 2X) , IX, Z2) 

WRITE (6 , 1000) DAY , HOUR, NBULK, DEN, NTHERM, BETA 
1000  FORMAT ( 1H  , 13, IX, 15 , 4F6. 1) 

WRITE  (9, 1025)  DAY,  HOUR,  DN 
1025  FORMAT (14, 14, IX, 33Z2) 

GO  TO  7 

185  WRITE (6, 1080) BLIM, ALIM 
1C80  FORMAT ('  BLIM  LE  ALIM  ',2E12.3> 

170  SUMLO=0.0 
180  SUMHI=0.0 
GO  TO  140 
7777  WRITE (6, 1050) 

1050  FORMAT (12HONORMAL  EXIT) 

STOP 

END 


I 


I 


1 CO  TO  189 

2 RETURN 

4 PRINT  K3CJMV7<J),R8<J>,N1<J>|H8<J,IC0) 

5 RETURN 

8 FRINT  * input  kappa 

9 INPUT  K3< J) 

10  RETURN 

12  PRINT  ht  factor  : 

13  INPUT  H3<J,K6> 

14  RETURN 


•1 


•I 


16  PRINT  •_  LIST  DATA  FROM  FILE  • *1 

< *»  1IA 


input  new  scale  : *1 


■l 


•I 


17  INPUT  NO 

18  FIND  M3 

19  GO  TO  83 
28  K8=l 
21  GO  TO  2188 

24  D2‘l 

25  RETURN 

28  PRINT  •_ 

29  INPUT  S 

30  .RETURN 

32  Le=2 

33  GO  TO  77 

36  L8=9 

37  RETURN 
48  D7-9 
41  GO  TO  189 

44  PRINT  ■_  input  den 

45  INPUT  NKJ) 

46  RETURN 

48  PRINT  •_  input  w.  i 

49  INPUT  KO<J> 

58  RETURN 

52  PRINT  •_  inpat  bulk  vz  : *| 
33  INPUT  U7(J) 

54  RETURN 
56  GO  TO  1810 

63  KG=2 

61  GC  TO  2178 

64  DS=0 

65  RETURN 

66  GO  TO  430 

68  se«i 

69  HOME 

70  PRINT  *nchan(0)a*| 

71  INPUT  Z8 

72  GO  TO  540 

73  GO  TO  1343 

76  CO  TO  1540 

77  07*1 

79  GO  TO  480 

80  07=1 

81  GO  TO  488 

82  GO  TO  ICO 

85  INPUT  (?33:K8,N9,H,H6,M7 

86  PAGE 

87  PRINT  HI* 

88  PRINT  N6 

89  PRINT  H7 

90  FOR  I«1  TO  24 

91  INPUT  §33: H 

92  PRINT  • . 

93  NEXT  I 

94  FOR  I«1  TO  11 

95  INPUT  033: H 

96  PRINT  M 

97  NEXT  I 

98  CO  TO  98 

188  PRINT  *.  INPUT  SPECTRUH  HOUR 
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lie  INPUT  H 
128  j«e 

130  J«JM 

140  IF  L2< J)«N  THEH  520 

150  IF  J<N7  THEN  133 

166  PRINT  SPECTRUM  NOT  FOUND' 

170  RETURN 

180  DIN  F(24t 15) » U7( 15) » N3< 15> » K8< 15) »H8<15t 2) » N1(15>»A0< 15) 
190  DIN  U(24),L1(16),L2(16> 

200  H0«1 

2ie  FOR  K«1  TO  24 
220  READ  U<K> 

230  NEXT  K 

240  DATA  0.167,0.189,0.209,0.227,0.243,8.259.0.273,0.287 

250  DATA  0.3,0.313,0.33,0.359,0.39,0.419,0.447,0.473 

268  DATA  0.497,0.52,0.543,0.565,0.586,0.666,0.626,0.645 

270  FOR  1*1  TO  15 

260  V7<I)=0 

290  NEXT  I 

3C0  D2=8 

310  N7=0 

320  S*1 

330  $0*8 

340  10*1  • 

350  DELETE  1,1 
366  PAGE 

370  UIHDOW  0,25,0,33 
380  UIEUPORT  30,133,20,189 
390  PAGE 

400  PRINT  'C  READY  ..  DATA  COMPARISON  PROGRAM...  G* 

418  RETURN 
423  K*0 
430  K*K+1 

440  IF  K>N7  THEN  478  , 

450  IF  K>K6  THEN  47A 


460  IF  N2(K)OL2(J>  THEN  430 
470  GO  TO  633 

480  PRINT  INPUT  SPECTRUM  INDEX  1 "I 
493  INPUT  J 
500  PAGE 

519  D7“l 

520  PAGE 

530  GO  TO  530 

540  M2*<1403.7TU(2S>-B>t2/At2 
550  Fl  = (<l+l/C)/<l+W2/C»tC*ietD 
560  U9=F<28, J>-F1 
570  D?=l 

580  UIEk'PORT  30,133,20,183 
590  IF  L0=0  THEN  620 
600  AXIS  2.137,6 
610  GO  TO  630 
620  AXIS  2.137,0 
630  Hl*l 

640  FOR  1*1  TO  24 
650  IF  L8=2  THEH  680 
660  X*30tUCI> 

670  GO  TO  690 
680  X*2C*LGT<3OWI»-10 
690  IF  £0=0  THEN  740 
7G0  F1=F(I, J)-F1-U9 
710  IF  FI >0  THEN  758 
720  Y*8 
730  GO  TO  859 
740  F1*F(I, J) 

750  Y=4.5tSTFl 

760  IF  L0=0  THEN  780 

770  Y*6*LGT(F1>+21 

788  IF  D7*l  THEH  S40 

750  MOUE  X, Y 

809  RMOUE  >0.13,-0.45 
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810  IF  Y<0  THEN  870 

829  PRINT 

830  GO  TO  870 

840  IF  I>1  THEN  880 
850  MOVE  X,Y 
860  DRAM  X* Y 
£70  NEXT  I 

III  IF  L0O1  THEM  1390 

9C0  IF  06=0  TKEH  13S0 
910  FOR  1*1  TO  24 
920  R2=I*11 
930  R3*1>11 
940  R1X1+R2+2JR3 

950  MOVE  30  JV< I ) , 6SLGT<0. 0069396166/R1 )+21 

960  RilOVE  -0.13,-8.45 

970  PRINT 

980  NEXT  I 

990  GO  TO  1390 

1030  RETURN 

1018  PRINT  RETRIEVE  DATA  FROM  FILE  • *J 

1020  J=1 
1033  INPUT  N9 
1643  FIND  N0 

1050  PRINT  * .START  WITH  SPECTRUM  DAY  HOUR:  "I 
1068  INPUT  M0, Ml 

1070  PRINT  • .LAST  SPECTRUM  AT  DAY  HOUR:  *1 

ices  INPUT  01,02 

1690  INPUT  033:19, L7 

1103  L9=L9*lSeS 

1110  INPUT  <?33:L1<J>,L2(J> 

1120  IF  LKJXrlO  THEN  1150 
1133  IF  L2<  JXM1  THEN  1150 
1149  GO  TO  1260 
1153  FOR  1*1  TO  33 
1160  INPUT  833! K7 
1170  NEXT  I 
1180  GO  TO  1119 
1193  INPUT  033: LI C J> , L2< J> 

12C3  INPUT  833: K7 

1210  FOR  1*1  TO  24 

1220  INPUT  e33:K7 

1239  R2=I*11 

1243  R3=I>11 

1250  Rl=l+R2+2tR3 

1260  F< I , J) =68. 861250. 964599tK7/Rl 

1270  NEXT  I 

1280  FCR  1=1  TO  8 

1290  INPUT  033: K7 

1309  NEXT  I 

1310  J=J+1 

1320  IF  J>16  THEN  1350 
1330  IF  LKJ-1X0I  THEN  1199 
1348  IF  L2U-1X02  THEN  1190 
1350  N7=J-1 
1360  C5=CNR<65+L7> 

1370  PRINT  DATA  RETRIEVED  GFROM  FILE  • "}N0 
1380  RETURN 
1390  HONE 

1400  PRINT  \._",L9 
1416  PR1HT  * Solrod  11"JQ* 

1420  PRINT  L1<J)5L2(J) 

1430  PRINT  " INDEX  * "IJ 
1440  IF  LQ=0  THEN  1490 
1450  IF  L0=1  THEN  1480 
1460  PRIHT  •_  Log-log  plot* 

1478  GO  TO  1490 
1480  PRIHT  ",  ScRi-lo3  Plot* 

1490  IF  S=1  THEN  1510 
1568  PRINT  • scale  ■ "IS 
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1510  HOME 

1529  L0-1 

1530  RETURN 
1548  PAGE 
1550  L=-l 
1560  K=1 

1570  PCS  XG=3  TO  99  STEP  33 
1560  FOR  YO=0  TO  75  STEP  25 
1550  UIENFCRT  X3,XS+30,75-Y0,95-Y0 
16C0  L=L*l 

1610  IF  L>0  THEN  1728 
1628  PRINT  •„  Solrad  ll'jQI 
1630  PRINT  •.  " J19 

1640  PRINT  •_  MLKl) 

1650  IF  L0=0  THEN  16e0 
1660  PRINT  •_  Seeii-loa  plot- 
1670  GO  TO  23e3 
1683  IF  S»1  THEN  2888 
1690  PRINT  ■ scale  ■ *» 

17C0  PRINT  S 
1710  GO  TO  2C30 
1720  IF  L>H7  THEN  2830 
1730  IF  L6=G  THEN  1760 
1740  AXIS  2.137,6 
1750  GO  TO  1773 
1760  AXIS  2.137,0 
1770  HOUE  15,25 
1780  FRINT  L2(L> 

1790  Hi  = l 

1888  FOR  1*1  TO  24 
1810  IF  L0=2  THEN  1840 
ie20  X=33tU(I> 

1833  GO  TO  1350 

1840  X*20$LGT<30WI>>-1O 

1850  F1*F<I.L> 

1860  Y=4.5iSTrl 
1878  IF  LC=3  THEN  1893 
18S3  Y=6SLGT(Fl>+21 
1659  IF  D7=l  THEN  1940 
19CQ  IF  Y<0  THEN  1970 
1918  ROUE  X,  Y 
1928  RERAN  8,0 
1933  GO  TO  1970 
1940  IF  I>1  THEN  1968 
1950  HOUE  X, Y 
I960  1 X,Y 

1970  NEXT  I 
1983  K=X*1 

1990  IF  L0O1  THEN  2080 
2030  IF  D8=3  THEN  2083 
2018  FOR  1=1  TO  24 
2020  R2=I=11 
2030  R3=I >11 
2040  Rl=l+R2+2tR3 

2e59  HOUE  30?U<I>,6*LGT(0.0969396166/R1H21 

2060  RDRAW  0,0 

2070  NEXT  I 

2080  NEXT  Y0 

2390  NEXT  X0 

2160  HOME 

2110  FRINT  'G-J 

2120  RETURN 

2130  INPUT  28 

2140  U2=<1403.7<UCZ8)-B)t2/AT2 
2150  Fl  = (<l  + l/C)-'(l+M2/C))tC»10'fD 
2160  U9=r (28, J)-F1 
2170  C=K8< J) 

2180  B=U7(J) 

2198  A=W0( J> 

2208  D=H8(J,K0> 
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2210  H*N1( J) 

2229  L=8 

2230  Ul=(B-3.1*A>/<6.794 
2240  U2=(B+3. 1*A)/4S.  794 
2250  IF  K0=1  THEN  2310 
2260  Ul-Ul 

2270  U2CU2 

2280  IF  Ut>3  THEN  2319 
2290  Ul=5 
2390  U2=23 

2310  FOR  X*U1  TO  U2  STEP  0.2 
2320  N2=(46.7S4JX-BK2/At2 
2330  IF  KC=2  THEN  2360 
2340  Ft=Ht£XP(-H2)*10tD 
2350  GO  TO  2120 

2360  Fl=<(l+l/C)/<l*H2/C)>tCll0tD 

2370  IF  LG=1  THEN  2480 

2380  Y=4.5:tS*Fl 

2359  GO  TO  2410 

2488  Y=6$LGT<F1>*21 

2410  IF  L<=0  THEN  2440 

2420  DRAM  X,Y 

2438  GO  TO  2450 

2440  HOME  X,Y 

2450  L=F 1 

2460  NEXT  X 

2470  PRINT  “UJJJJJJJJJJJJ* 

2480  IF  K0=2  THEN  2510 
2420  PRINT  ■Maxwellian  fit* 

2500  GO  TO  2520 

2510  PRINT  “Kappa  dist  fit' 

252e  PRINT  '.Sulk  Uz  =•  -;U7<J> 
2530  PRINT  "Them  u.  = “}A 
2540  PRINT  “Den  = “iNKJll Vcc“ 
2553  IF  K0S1  THEN  2570 
2560  PRINT  “Kappa  ■ "JC 
2578  RETURN 


1 GO  TO  200 

2 RETURN 

4 LO=0 

5 GO  TO  390 

8 PRINT  *>0  ALPHA  EQUALS  *| 


10  INPUT 

11  GO  TO 


A 
390 


390 

•0(X,Y> 


■SCALE  EQUALS 
S 


12  PRINT  "NEH  BETA  EQUALS  "1 

14  INPUT  B 

15  GO  TO  390 

16  PRINT  *_  LIST  DATA  FROM  FILE  • *J 
1?  INPUT  H0 

18  FIND  N0 

19  GO  TO  85 

20  PAGE 

21  LIST 

22  RETURN 

23  GO  TO  390 

24  IF  H0=0  THEN  1240 

25  GO  TO  1270 
32  GO  TO  180 

36  L0=9 

37  GO  TO  41 

40  L0=1 

41  PRINT  ‘INPUT  SPECTRUM  INDEX:  *J 

42  INPUT  J 

43  GO  TO  1190 

44  L0=1 

45  GO  TO 

48  PRINT 

49  INPUT  X9.Y9 

51  GO  TO  398 

52  PRINT 

53  INPUT 

54  GO  TO  390 
56  GO  TO  1600 
60  GO  TO  1370 

64  H0*1 

65  GO  TO  360 
72  GO  TO  100 
76  GO  TO  360 
88  C0=1 
81  GO  TO  410 

85  INPUT  B33:H3»N9»M»M6»N7 

86  PAGE 

87  PRINT  tV 

88  PRINT  N6 

89  PRINT  M7 

90  FOR  1=1  TO  24 

91  INPUT  033 :M 

92  PRINT  ■ *»H 

93  NEXT  I 

94  FOR  1=1  TO  11 

95  INPUT  033 :M 

96  PRIHT  M 

97  NEXT  I 

98  GO  TO  90 

100  PRINT  INPUT  SPECTRUM  HOUR 

110  INPUT  M 
120  J=0 
130  J=J+1 

140  IF  L2< J)=I1  THEN  180 
150  IF  J<N7  THEN  130 
168  PRINT  *_  SPECTRUM  NOT  FOUND* 

178  RETURN 
180  L0=0 
198  GO  TO  1190 
200  DIM  F <24 * 20 ) 

210  DIM  C*<2)iA$<I).B»Cl) 
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220  Din  U(24),T1(30)|L1(30),L2(30) 

230  H0=1 
240  N7=0 
250  DELETE  1,1 
260  FOR  K»1  TO  24 
278  READ  U(tO 
200  NEXT  K 

290  DATA  0.167,0.169,0.209,0.227,0.243,0.259,0.273,0.287 

300  DATA  0.3,0.313,0.33,0.359,0.39-0.419,0.447,0.473 

310  DATA  0.497,0.52,0.543,0.565,0.586,0.606,0.626,0.645 

320  READ  A,B,S,T2,E1,E2,X9,Y9 

330  DATA  30,45,2,0,2,2,0,0 

340  PAGE 

350  SET  DEGREES 

368  PAGE 

370  PRINT  'G  READY  _ PROGRAM  DOES  NOT  INCLUDE.  KEYBOARD  DN  INPUT. 
3ee  RETURN 
390  PAGE 
400  C0=0 

410  WINDOW  -12. 5-X9 , 52 . 5-X9 , -20- Y9 , 30-Y9 
420  VIEWPORT  0,130,0,100 
430  H0=0 
440  H1=0 

450  X0=20*S*COS<B) 

460  Y0=20tS:SIH<B)*SIN<A) 

470  Yl=12*StC0S<A> 

4e0  X2=20tS*SIN(B> 

490  Y2=-20*S*C0S<B)*SIN<A> 

580  I10UE  0,0 
510  DRAW  X0,Y0 
520  PRINT  • t* 

530  MOVE  0,0 

540  DRAM  0 , Y 1 

550  IF  LOal  THEN  580 

560  PRIHT  ’ ('  " 

570  GO  TO  598 

589  PRINT  " log  f* 

590  MOVE  0,0 
600  DRAW  X2»Y2 
610  PRINT  * v" 

620  FOR  K=1  TO  24  STEP  23 
630  X3=X2*V<K> 

640  Y3»Y2*V<K> 

650  MOVE  X3,Y3 
660  RDRAM  XO,Y0 
670  NEXT  K 

680  IF  C0*1  THEN  930 
690  T2*0 
700  T1(N7)=0 

710  FOR  J=1  TO  N7  . . . 

720  IF  J=1  THEN  760 
730  T2*T2+T1 (J-l ) 

740  T*J+T2-1 
750  GO  TO  770 
760  T=0 

770  FOR  1=1  TO  24 
780  V1=20*'KI> 

790  D=S0R(Ult2+Tt2> 

800  C=B*AT1!<T/U1) 

818  x*s:d*sin<c> 

820  IF  LO=0  THEN  860 
838  Y8=3*LGT<.F<I,J>'*8 
840  Y«-SmC0S<C)*SIN<A>*S*C0S<A)*Y8 
850  GO  TO  870 

660  Y«-StD*COS(C)»SlN(A>+StCOS(A)»F(l,J) 

870  IF  I>1  THEN  890 
880  MOVE  X, Y 

690  DRAW  X,  Y .'  • 

980  NEXT  I 
910  NEXT  J 
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920  CO  TO  2050 
930  FOR  1*1  70  24 
940  V1=20*U(I) 

950  V2=Vlt2 

960  T3=0 

970  T 1 (H7>*0 

980  Y8=-l 

990  FOR  J*1  TO  N7 

1000  IF  J=1  THEN  1040 

1010  T3=T3*Tl<J-l> 

1020  T=J+T3-1 
1030  GO  TO  1050 
1040  T=0 

1050  D=SGR<U2*Tt2> 

1060  C=B+ATN(T/V1> 

1070  X=S*D*SIH<C> 

1080  IF  L0=0  THEM  1120 
1090  Y6=3*LGT<F(I,J)>+8 
1100  Y=-S*DtCOS(C) JSIN(A>+S*C0S(A)>Y8 
1110  GO  TO  1139 

1120  Y=-S*D*COS<C>t$IN<A>+StCOS<A>*F<I,J> 

1130  IF  J>1  THEH  1150 

1140  MOVIE  X,Y 

1150  PRAM  X,Y 

1160  NEXT  J 

1178  HEXT  I 

lieO  GO  TO  2050 

1190  WINDOW  0,25,8,30 

1200  VIEWPORT  30,130,20,100 

1210  PAGE 

1220  AXIS 

1230  Hl*l 

1240  FOR  1*1  TO  24 
1250  X=20ttKl> 

1260  IF  L0*0  THEN  1290 

1270  Y=S*E1*(L0G(F<I, J))*E2) 

1280  GO  TO  1330 
1290  Y*S1F(I,J> 

1300  IF  I>1  THEN  1320 
1310  MOVE  X, Y 
1320  DRAW  X,  Y 
1330  HEXT  I 
1340  HOME 
1350  GO  TO  2050 
1360  RETURN 
1370  K?=0 

1380  PRINT  •_  STORE  DATA  ON  FILE  • 'I 

1390  INPUT  N0 

1400  FIND  H0 

1410  PRINT  033:L9 

1420  PRINT  033:17 

1430  FOR  K=i  TO  H? 

1440  PRINT  033:L1<IO 
1450  PRINT  033:L2OO 
1460  PRINT  033*. K7 
1470  FOR  1=1  TO  24 
14e8  R2=I=1 1 
1490  R3=I >1 1 
1500  R1=1+R2+2*R3 

1510  K0=1  17. 09429-27. ?44786*L0G(R1*F<I,K>> 

1520  PRINT  033: K0 

1530  HEXT  I 

1540  FOR  H«l  TO  8 

1550  PRINT  033: K7 

1560  HEXT  M 

1570  HEXT  K 

1580  PRINT  ’DATA  STORED  IN  FILE  • *1118 
1590  RETURN 

1600  PRINT  RETRIEVE  DATA  FROM  FILE  » *1 
1610  J*1 
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1629  INPUT  N9 

1630  FIND  N0 

1640  PRINT  * .START  WITH  SPECTRUM  DAY  HOUR:  *| 
1650  INPUT  HO, Ml 

1660  PRINT  ' .LAST  SPECTRUM  AT  DAY  HOUR:  "I 

1670  INPUT  01,02 

1680  INPUT  P33:L9,L7 

1690  L9=L9+1 900 

1780  INPUT  P33:L1<J>,L2<J> 

1710  IF  LKJXNO  THEN  1750 
1720  IF  L2CJXM1  THEN  1750 
1733  L8=L2< J) 

1740  GO  TO  1790 
1750  FOR  1=1  TO  33 
1760  INPUT  P33JK7 
1770  NEXT  I 

1780  INPUT  033:L1 (J),L2< J) 

1790  INPUT  033: K? 

1808  FOR  1=1  TO  24 
1810  INPUT  P33:K7 
1820  R2«I*ll 
1830  R3=I >1 1 
1840  R1=1+R2+2*R3 

1850  F(I,J)=68. 0612*0. 964599tK7/Rt 

I860  NEXT  I 

1870  FOR  1=1  TO  8 

1880  INPUT  P33JK7 

1890  NEXT  I 

1908  IF  J=1  THEN  1990 

1910  L6=L8 

1920  L8=L2( J) 

1938  L3=L8-L6 

1949  IF  L3>0  THEN  I960 

1959  L3=L3+2400 

1960  IF  L3<40  THEH  1980 
1970  L3»L3-40 

1980  Tl(J-l)=L3/2-l 
1999  J*J+1 

2800  IF  LKJ-1X01  THEN  1780 
2010  IF  L2CJ-1X02  THEN  1780 
2020  N7-J-1 

2030  PRINT  DATA  RETRIEUED  GFROM  FILE  t "IHO 
2040  RETURH 
2050  HOME 

2060  PRINT  " ",L9 

2078  IF  L7=l  THEN  2100 
2080  PRINT  • Sol  rod  UA* 

2090  GO  TO  2110 

2100  PRINT  ’ Solrad  11B" 

2118  IF  Hl=l  THEN  2150 
2120  PRINT  LK1>;L2(1) 

2130  PRINT  LKN7)!L2(H7> 

2148  GO  TO  2198 
2150  PRINT  L1<JXL2<J> 

2160  PRINT  • INDEX  * "JJ 
2170  IF  LO=0  THEN  2190 
2180  PRINT  ■_  Scni-loa  plot" 

2190  HOME 
2200  RETURN 
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